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f{R) gravity is one of the simplest theories of modified gravity to explain the accelerated cosmic 
expansion. Although it is usually assumed that the quasi-Newtonian approach (a combination of 
the quasi-static approximation and sub-Hubble limit) for cosmic perturbations is good enough to 
describe the evolution of large scale structure in f{R) models, some studies have suggested that 
this method is not valid for all f{R) models. Here, we show that in the matter-dominated era, the 
pressure and shear equations alone, which can be recast into four first-order equations to solve for 
cosmological perturbations exactly, are sufficient to solve for the Newtonian potential, and the 
curvature potential, $. Based on these two equations, we are able to clarify how the exact linear 
perturbations fit into different limits. We find that the Compton length controls the quasi-static 
behaviours in f{R) gravity. In addition, regardless the validity of quasi-static approximation, a 
strong version of the sub-Hubble limit alone is sufhcient to reduce the exact linear perturbations in 
any viable f{R) gravity to second order. Our findings disagree with some previous studies where 
we find little difference between our exact and quasi-Newtonian solutions even up to fc = lOc~^'Ho. 

PACS numbers: 98.80.Jk 


I. INTRODUCTION 

Among all the possibilities to explain the observed accelerated expansion of the Universe with a modified theory of 
Einstein’s General Relativity (GR), f{R) gravity (also dubbed fourth-order gravity) is the simplest one EH- In this 
class of theories, a new function of the Ricci scalar, f{R), is included in the Einstein-Hilbert action. In this way, f{R) 
models form a class of higher derivative gravity theories. This is a natural extension of Einstein’s General Relativity 
because there is no prior reason to exclude these higher order terms from the Lagrangian density Indeed, the 

higher order terms of the Ricci scalar, R, does appear in low energy effective Lagrangians in string theory and other 
candidate theories of quantum gravity @ . 

Phenomenally speaking, there are viable f{R) models that can not only yield a consistent and realistic cosmology [l3, 
[H , but also pass the solar system test 0- However, although some /(i?) models can produce a background evolution 
that is identical to the standard cosmology (AGDM), this same set of f{R) models will differ from AGDM with regard 
to other phenomena, such as weak leasing, cluster abundance, cosmic microwave background (GMB), and baryon 
acoustic oscillations (BAG) [l^ - [l5l| . To understand these phenomena on large scales, cosmological perturbation 
theory is a significant building block. According to this theory, the original fluctuations were amplified beyond 
the scale of the horizon at the the end of inflationary epoch, and only after the horizon has grown to the size of 
the fluctuations, various structures we recognize today start to grow (see e.g. 0). Moreover, we also learn from 
cosmological perturbation theory that there are three types of perturbation: tensor, vector, and scalar modes Ea. Of 
these, the tensor perturbation will cause gravitational wave, and the vector perturbation will generate vorticity that 
decays with time and becomes entirely negligible. This leaves the scalar perturbation as the only source that could 
contribute to the growth of cosmic structures. 

Many papers have been devoted to develop the theory of cosmological linear perturbations within the framework 
of /(i?) gravity EB“I13' However, many of them si mply formulate the linear perturbations in f(R) models by taking 
the quasi-static assumption for granted E^, [l^ [2^ - 1?^ . Even if several works do provide us with ways to solve for 
linear perturbations exactly, it is not always clear under what circumstance the exact linear perturbations could 
be approximated by the quasi-static solutions well [20l - [^ [^ . Indeed, some studies suggest that the quasi-static 
approximation will break down outside the sound horizon of modified gravity [s^, or even, for certain f{R) models, on 
sub-Hubble scales [2^; others conclude that the quasi-static approximation is valid for the most practical /(i?) models 
either on the nonlinear [3l|, or near-Hubble scales [s^. To clarify the reason behind these seemingly contradictory 
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results, analyzes the quasi-static approximation in f{R) models within the Einstein frame, and find that the 
fast/slow-rolling behaviours of the cosmic background will determine the validity of the quasi-static approximation. 

In this paper, we clarify when the quasi-static approximation will break down within the framework of metric f{R) 
gravity alone. In other words, unlike [s^, our analysis is purely based on linear perturbations of f{R) gravity in the 
Jordan frame, and will not invoke the ideas of a fast/slow-rolling background as the explanation. For this purpose, we 
firstly formulate the exact linear perturbations in a simple and clarifying way, and then investigate how these exact 
linear perturbations fit into different limits analytically as well as numerically. 

The structure of the paper is organized as below. In Sec. [Ill we briefly outline the basic formalism of f{R) gravity, 
and its application to the cosmic background evolution. We then show how to reduce the modified Einstein equations 
and equations of conservations for /(i?) models into the pressure and shear equations from first principles, and prove 
analytically why these two equations are sufficient to solve for the Newtonian potential and spatial curvature. In 
Sec. [nil we describe how linear perturbations in f{R) gravity fit into various limits, and discuss the role of the sub- 
Hubble and the quasi-Newtonian assumptions in taking the limits. In Sec. IlYl we apply the numerical method to 
solve for the Newtonian potentials and spatial curvatures exactly from the pressure and shear equations, and compare 
the solutions with those obtained in various limits. We also use numerical solutions from our newly derived exact 
equations to rebut the concerns over the quasi-static approximations. Finally, we outline a brief conclusion in Sec. El 


II. FORMALISM OF f{R) GRAVITY 
A. Background evolution in f{R) theories 


In /(i?) gravity, the Ricci Scalar, R, in the Einstein-Hilbert action is generalized to a function of R, and the action 
can be written as 


S = 



R + f{R) 

16ttG 



( 1 ) 


where is the Lagrangian for matter. According to this action, the effect of f{R) models can be understood as an 
additive correction to Einstein’s standard gravity (GR; we will use GR and standard gravity interchangeably), so it is 
easy to understand how f{R) gravity deviates from AGDM. This action also implies that there is a scale to distinguish 
the dominated terms between Einstein’s gravity and f{R) theory. In this notion of f{R) gravity, the modified field 
equations will be 

Gfiv + IrR^lv — (tIr = ( 2 ) 

where fn = df{R)/dR. 

If we only consider a flat, isotropic, and homogeneous Universe, it is still legitimate to apply the Robertson-Walker 
(RW) metric to f{R) cosmology. Under this metric, the trace of the Ricci tensor alone will yield the equality, 

R = + 6a-^n' = (3) 


which is purely geometric and independent of the choice of theory of gravity (hereafter, we define “ ' ” for d/da with 
R = a/a = dlna/dry, where a dot is denoted as d/drj and ry is conformal time). This equality tells us that as long 
as the RW metric is assumed, there will be one simple relation between the Hubble rate, H(a), and the Ricci scalar, 
R{a). Eq. ([3]) also tells us that the Ricci scalar is spatially invariant in an isotropic and homogeneous Universe; thus, 
only the temporal components of and afR in Eq. ([3]) will be nonzero. Indeed, like the Friedmann 

equation in Einstein’s standard theory of gravity, the 0 — 0 component of Eq. alone is sufficient to yield a modified 
Friedmann equation. 

Accordingly, after combining all the temporal components in the field equations, Eq. (U), we will obtain the modified 
Friedmann equation in /(i?) theory, 

Ite - afRUn' + \fc? + U^fRRaR! = ^a^p. (4) 

6 3 

From this modified Friedmann equation and Eq. ®, the Hubble function can then be expressed as a function of the 
Ricci scalar, 


SttGo^p/?, - a^{f - /flR)/6 

1 + /fl + Jrro-R' 


( 5 ) 
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where p = ‘i'Hfflm !in the matter-dominated Universe, and Hq is equal to the current value of the Hubble rate 
in ACDM. For clarity, unless specifically noted, we choose units to have c = Hq = 1 hereafter. 

In order to tell how the background of the Universe in f{R) gravity evolves with time, we need to solve for the 
Ricci scalar (or Hubble rate) as a function of time (or scale factor) explicitly. For this purpose, we find the most 
intuitive way is to recast Eq. (|4]) into a second order equation for the Ricci scalar. After differentiating Eq. ([5]) with 
respect to a, and combining it with Eq. (|4]) and Eq. m to eliminate 'H'H', we obtain a new nonlinear second order 
differential equation; 


R" + 


( f RRR I 


aR 



K^fRRR^ J 


R “h 


SirGp 

JrrR^ 


4(1 + Jr) 
o^Jrr 


( 6 ) 


where is given by Eq. Although Eq. ([6|) is mathematically e^ivalent to other equations written down in 
previous work for the cosmic expansion in f{R) theory (see e.g. [1, ® has the advantage that it shows 

the evolution of the Ricci scalar as an explicit solution of a second order equation. Like all second order differential 
equations, a negative coefficient of the second term in Eq. (|6]) will lead to a positive feedback system, which is highly 
unstable if the sign of this coefficient never changes. The role of this term becomes even clearer when the stability of 
f{R) theory is considered [s^. 


B. Cosmological perturbations in conformal Newtonian gauge 

For any covariant linear perturbation around a fiat, isotropic and homogeneous Universe (see e.g. Hil) , we can 
decompose the perturbation into scalar, vector, and tensor modes, with a harmonic expansion on a 3D sphere (for a 
4D universe). Because only scalar perturbations will contribute to the growth of structure, we do not consider vector 
and tensor perturbations in this paper. In addition, we will only focus our analysis on the conformal Newtonian 
gauge, which will yield the Newtonian-like equations, and is popular in the literature of structure formation and weak 
leasing [s^. For linear perturbations in a generic gauge in f{R) gravity, we refer the readers to M- 
Following the notation of [1^ , we write down the RW metric in the Newtonian gauge as 

^ 2^)dr]^ + 0^(1 - 2d>)da;^ (7) 

where 4' is called the Newtonian potential and <I> is the spatial curvature. We then apply this metric to the conservation 
equations, and obtain 

S+{l + w)i0-S^) + 3n{cl-w)6 = O, (8) 

e + n(i-3w+—^^e-^^^6-k'^'ii = o. (9) 

\ 1 + wj l-fw 

Here, 6 is defined as the density perturbation, 9 as the amplitude of spatial velocity, the speed of sound is = 6P/6p, 
and P/p = w is the equation of state. Similarly, we can apply the metric, Eq. ([71), to the perturbed modified Einstein 
equation. Accordingly, we will derive the Poisson equation, 

(l + /fl)[-fc^($ + 4')-3'H($ + 4') + (3?7-6-H2)5-_3-H$] + /)^(_9-H^f+3-H$_3$) = 8TrGpa^S+ik‘^-3n)a+3na, (10) 

the pressure equation, 

(l + /fl)[^ + $ +3-H(4'-K$)-1-3774' +(77-f 277^)$] 

+ /ii(3774' - 77$ + 3$) + ^(3$ - $) = SnGa^SP + (277^ + 77 - ^k^^ ct + 77a - d, (11) 

the momentum equation, 

(1 + /«)[$ + $+ 77(4- + $)] + ^(24- - $) = SiTGpa'^e/k'^ +'Ha-d, (12) 

and the shear equation, 

(4- - $) + —^ x [-6(77 + 77^)4- - 3774' + ^^4' - 977$ - 3$ - 2fc2$], (13) 

1 + /fi a (1 + Jr) 

where 5P is the perturbation of pressure density. We also define 

d = 12'KGa^{p + P)a/fc^, 


(14) 
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where a is anisotropic stress. Unlike the their counterparts in GR, in f{R) theories the anisotropic stress not only 
exists in the shear equation, but also appears in Poisson equation, the pressure equation, and the momentum equation. 
However, when f{R) is close to zero, we can eliminate the anisotropic stress terms from the Poisson, the pressure, 
and the momentum equations by substituting a in the shear equation. This will make a only appear in the shear 
equation, and recover the GR regime. It is also worth noting that in these equations, 

Ir = /rrR, ( 15 ) 

Ir = IrrrR^ + JrrR, ( 16 ) 

where Jrr = d? f /dR? and Jrrr = d^f/dR^] hence Jr, Jrr, and Jrrr all affect the linear perturbations in f{R) 
gravity in an explicit way. 

Since so far we have only assumed a flat, isotropic and homogeneous Universe, the conservation equations, Eq. ([5])- 
and modified Einstein equation, Eq. (nnii-dTsi), can also be applied to the radiation-dominated era as well as the 
matter-dominated era. Although it is quite straightforward to derive Eq. (ITUl) - (fO]) . they do not seem to have appeared 
in the literature previously. 


C. Matter-dominated era 


Considering the success of the standard cosmological model (ACDM) on the cosmic microwave background 
(CMB) [ 33 ,[s^, we follow to take the assumption that f{R) models will recover GR before the CMB is formed. 
Accordingly, we will only analyze the linear perturbations in the matter-dominated era. 

In the matter-dominated era, radiation density is ignored compared to matter density, so that cj. = SP/6p ~ P/p = 
w = 0, and the conservation equations, Eq. ®-®, can be reduced into 

d + nS + k"^^-3n^-3^ = 0. (17) 


In addition, it is also usually assumed that the Universe is free of anisotropic stress in the matter-dominated era, 
so it is legitimate to reduce the modified Einstein equations, Eq. (linD-disi), further by setting cr = 0. IndeedUi we 
take (7 = 0, the modified Einstein equations, Eq. (nnii-dTsi), will be identical to the field equations derived by [2^ to 
solve for matter density perturbation. However, what differs from is that we do not attempt to recast all six 
equations (2 conservation equations and 4 fields equations) into one single differential equations of matter density 
contrast. As we are going to show, it will become a much easier task if we solve for 4* and first, and then put these 
two potentials back into the Poisson equation and the momentum equation to obtain evolutions of S and 0, because 
the only equations required to solve for both potentials are the pressure equation. 


(1 + /ii)[^ + ^ + 3-H(4' + ^)+ 3774' + in + 2772 )$] -b /(r( 37^4' - 7^$ -b 34') -b /fl(34' - $) = 0, 


(18) 


and the shear equation, 

(4' - $) = — x [-6(77 -b 7^^)4' - -b - 97^$ - 3l> - 2k'^^]. (19) 

a (1 + /fl) 

To show the redundancy of the conservation and fields equations in f{R) gravity, we firstly combine the conservation 
of momentum equation, Eq. ®, and the momentum equation, Eq. m, by eliminating 6 in both equations. We will 
thus obtain 


(1 + /fl)[^ + ^ + 37^(4- -b 4>) -b 77(4- + $) -b 277^(4- -b 4>)] -b /(r( 5774' - 774> -b 34-) -b /fl(24' - $) = (20) 

Then it is straightforward to prove that the difference between this equation and the pressure equation, Eq. dn, 
tells us no more than how the cosmic background evolves. In other words, in terms of the evolutions of 4* and 4>, 
information given from the conservation equation for momentum and the momentum equation is equivalent to that 
derived from the pressure equation alone. Because of this equivalence, we call Eq. (OUl) the “alternative pressure 
equation”. Similarly, by combining the conservation equation for the density perturbation, Eq. ® , Poisson equation, 
Eq. (fTHl) . and the shear equation, Eq. (IT^ . we can cancel out 5, and show that a combination of these three equations 
does not tell us more than Eq. (fTOl) (the shear equation). Hence, we conclude that the pressure equation, Eq. (fTSl) . 
and the shear equation, Eq. m, alone are sufficient to solve for 4> and 4^. The “alternative pressure equation”, 
Eq. dini), show explicitly how matter density will affect evolutions of both potentials. On the contrary, if we choose 
the pressure equation, Eq. (HID, over the “alternative pressure equation”, Eq. (EHl), then the role of matter density 
will become implicit; that is, the evolutions of both potentials will be determined purely by the Hubble rate, 77, and 
the Ricci scalar, R. 
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III. COMPTON LENGTH AND EFFECT OF f{R) GRAVITY 


In the previous section, we have shown that the coupled pressure and shear equations are all we need in order to 
solve for d/ and $ exactly. To understand how these equations fit into various limits of cosmological perturbations in 
f{R) theories, we define a scale length. 


= 


Irr 


^(1 + / h )’ 


( 21 ) 


which is easily recognizable in the shear equation, Eq. m- Based on this scale length, the elements of the shear 
equation can be classified into two groups: one is proportional to the scale dependent 1+k^Xj] the other is proportional 
to just Aj, where k~^ and A/ are in the unit, T-Lq/c. It is worth noting that although derived from a different set of 
equations. A/ looks very similar to the Compton wavelength in 0, which we compare with in Appendix |B1 or the 
lengthscale defined in 

Before we discuss how these factors will affect linear perturbations in f{R) gravity, we would like to clarify our 
definition of the so called quasi-static approximation, which is sometimes confusing in the existing literature (cf. 
[2i,[3l|3l,[33). According to our definition, the quasi-Newtonian approximation contain two parts: the sub-Hubble 
limit ( k ^ R ) and the quasi-static approximation (|A| < 'H|A| ), where X might be H, T or <I>. Not following the 
notations in [U, [H, , we call k ^ R as the sub-Hubble limit rather than the sub-horizon limit because 

R is not a horizon but the Hubble length, as well as because we would like to emphasize the difference between this 
lengthscale and the particle horizon of modified gravity discussed in [s^ . 

It shall be kept in mind that although R^ might become large at the early epoch, its effect must be balanced out 
by Xj, otherwise f{R) gravity will not match ACDM in the early Universe. By contrast, when the effects of f{R) 
gravity start to kick-in at the more recent epoch, R should have become compatible with the order of Rq. Based 
on this fact, we also introduce a weaker version of the sub-Hubble approximation, k 3> Ro- Any k that satisfies the 
sub-Hubble limit must also satisfy this weaker version of the sub-Hubble limit, but not vice versa. 


A. Reduced second-order differential equations 


Without the need to invoke the quasi-static approximation, when 1 -|- fc^A^ ^ can reduce the shear 

equation, Eq. CM, into a simple relation between T and $ by dropping terms proportional to A^. Here we introduce 
a wavenumber, 

= max{R^, RX/X, X/X}, (22) 


which is defined by the maximum value of expansion rate of Universe or rate of change of X. Again, X is defined as 
R, T or however, our discuss below will only focus on X as representing tj/ or $, because for any f{R) model that 
yields a similar cosmic acceleration to ACDM, R^ ~ R^ R. In terms of the scale length. A/, Eq. (1191) reduces to 

{l-^4k^Xj)^ = {l-^2k^X})^. (23) 

Note that this immediately gives us the gravitational slip between the potentials, This relation then can be put 

back into the alternative pressure equation, Eq (1^ . to yield the reduced second-order differential equations. 


^ -f {m + Ai- As) i> + (r + 2R'^+A 2 + = 0, 

$ + {3R + Ai+ A 3 ) i>+(R + 2R^ + A2- = 0, 

where we define 

. 3//j(l-f 4fc^Aj) 

' = (l + /fl)(2 + 6A:2Ap’ 

_ + Jr {mh + 3/fl) (1 + Ak^X})aHTTGp 

(1 +/fl)(2-f 6A:^Aj) (1-f/_ r)( 1-f 3fc^A^) (1-I-3fc^Ay)(l-I-/_ r) 

^ 4k%Xf 

^ (l + 2FA2)(l + 3fc2A2)’ 


(24) 

(25) 


(26) 

(27) 


( 28 ) 
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2P(4fc2A/A2 - snXfXf -X}- XfXf) 

(l + 4fc2A2)(l + 3fc2Ap 

2k'^{8k'^Xf^ kj — yUXfXf — kj — XfXf) dk"^ fiiXfXf 

(l + 2fc2A2)(l + 3fc2Ap (1 + /fl)(l + 2A:2A2)(1 + 3k^Xj)' 


(29) 

(30) 


In these equations, //{, fu, A/, and A/ can be derived from background evolution without solving Eq. (l24ll and 
Eq. (l25ll . In order to get a sense about the order of magnitude of these terms, we take the quasi-static approximation 
for R and A/, and obtain Jr ~ /rrH^R, Jr ^ {Jrrr -t fRR)'H^R, A/ ~ RXf, and A/ ^ 'H^A/, which we shall keep 
in mind are only approximately correct. In addition, for any viable f{R) model, the function f{R) can only grow 
steadily from an asymptotic constant with /_r < 0 and Jrr > 0, which implies that |/_r_r_r| must be less than /rr. 
Hence, when f{R) effects kick in, Jr and Jr are of the same order of magnitude as Xj. 

Based on the approximations above, it is reasonable to argue that the relationship between the exact but coupled 
second-order differential equations, Eq. (HU) and Eq. o, and the two reduced decoupled differential equations, 
Eq. (IMl) and Eq. (I25L depends on the competition between k^Xj: and 1 -|- k'^kj. Indeed, the condition, 1 -|- k'^kj ^ 
k\X^f, can be further broken down into A/ <C 1/fcjc or k ^ kx, so either of which will be sufficient for the validity 
of Eq. (IMl) and Eq. (OSl) . Henceforth, we shall conclude that this transition from the exact equations to the reduced 
second-order ones will happen as long as either f{R) effects become trivial (A/ <C l/fcjc) or the scale of the system 
is relatively small (fc ^ kx)- We call k ^ kx the strong version of sub-Hubble limit because this condition naturally 
contains the other two versions of sub-Hubble limit. In the limit. A/ <C i/kx, all terms proportional to A/ or less in 
Eq. (IMl) and Eq. (1^51) are far less than H, and f{R) effects will become significant only when kXf is nontrivial, that 
is, when A/ > 1/k. Similarly, in the limit k ^ kx, f{R) effects will dominate only when k > 1/A/. 

We should emphasize that although there are two separate conditions in our discussion for a given f{R) model, the 
evolution of 'k and might satisfy A/ 1/kx at the earlier time, and then satisfy k ^ kx only recently. If this is 
the case, and if time derivatives of d' and $ are negligible at later epoch, Eq. (IMl) and Eq. (OSl) will only require the 
weak version of the sub-Hubble limit, k ^ Hq, for their validity, because during this period, H. is compatible with 
Ho- 

So far we have overlooked what kx really represents by simply assuming that we can always find a k that is large 
enough (or A/ that is small enough) to guarantee k ^ kx (or Xf 1/kx)- In fact, we may or may not find a 
case where the ratio between the time derivatives of potentials and potentials themselves, which we denote X/X and 
X/X, are so large that we will never find a reasonable k or Xf for these approximations. To check this, we go back 
to analyze Eq. (IMl) and Eq. (OSl) themselves. If these equations yield highly unstable solutions (that is, 4' or that 
increases/decreases significantly within a short period of time), our assumptions might be flaw; otherwise, they are 
fairly good. It turns out that for models that match ACDM at early time, a positive coefficient of the first derivatives 
in Eq. (HI and Eq. (1^51) will always keep the solutions stable, which naturally includes the cases when A/ <C 1. 
Indeed, when A/ <C 1, the potentials will hardly grow faster than In a; thus, the quasi-static approximation will never 
break down. Even if |//i| A^ is compatible with TL, and makes the coefficient of the first derivatives in Eq. (IMl) and 
Eq. (ESI) negative (remember Jr < 0), while k is large enough, ^ HX/X or X/X will hardly break down because 
the potentials must grow as fast as to make X/X ~ k^/H. Perhaps the only way to achieve this is when 

Xf H, which clearly does not fit any viable f{R) model. 


B. Quasi-Newtonian approximation for matter density contrast 

Taking the strong version of sub-Hubble approximation, k ^ kx, alone, we can drop all terms that are not 
proportional to fc^, so the equation of conservation of energy-momentum, Eq. (II3, will be reduced to 

6 + n6 + k‘^^ = 0. (31) 

Similarly, the Poisson equation, Eq. (HOD, and the shear equation, Eq. dUD, will lead to a simple relation between 
Newtonian potential and the matter density perturbation (iSl Il9j|. 

2 -(- Qk^X^ 

Gfc2^ = 87rGpo^ and G =-a{lJ r)-^-^-^^^, (32) 


where po is a constant in the matter-dominated epoch. It shall be noticed that according to Eq. (15^ . 4' basically 
grows with a factor of (5/a; thus, even if 4^ behaves like a constant, time derivatives of 5 could still be large. 
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After putting Eq. (1321) back into Eq. dSU, we will obtain a second order differential equation for S, which often 
appears in the f(-R) literature because of the quasi-Newtonian approximation [H El, m m, m, [13. However, as 
we have just shown. Eg. (1311) and (1321) can be derived from the exact Poisson, shear and conservation equations, even 
if we take a more broad approximation, k ^ kx, which includes the quasi-Newtonian approximation as well as the 
limit, k ^ kx In other words, as long as k is large enough, Eg. dSTI) and (l32|) will always be correct, no matter 
if the quasi-static approximation breaks down or not. 

In order to compare Eq. d^Tl) and (1321) with the reduced second-order differential equations in the last section, after 
some manipulation, we rewrite Eq. (1311) in form of a second order differential equation of 'k. 


^ + (3n + Ai- A 3 ) i’+(n + 2-H^ + A2 + h^) = o. 


(33) 


where A 3 is defined in Eq. (|2^ . in Eq. (|29)) . and 


SUfR + fR il+4eX})aHTTGp 4fc2/flA/A/ 

1 + fR {l + 3k^X}){l + fR) (l + /fl)(l + 4fc2Ap(l + 3fc2A2)- 


(34) 

(35) 


A comparison between Ai, A 2 , and Ai, A 2 shows that when kXf ^ 1 ot /r 1, which also implies A/ <C 1, Eq. (l24ll 
and Eq. (IM)) will become identical. This is interesting because A/ <C 1 is included in the approximation. A/ <C Ijkx, 
which is one of the two conditions that will reduce Eq. (IMl) . This means that when A/ is small enough, Eq. (1551) . like 
Eq. (l24l) . will still be applicable even at near-Hubble or super-Hubble scales. This is in agreement with slow-rolling 
solutions in [s^ • On the contrary, when A/ is not so small, the condition kXf ^ 1 will be consistent with the strong 
version of sub-Hubble limit, k ^ kx] thus, not surprisingly, Eq. (1241) and Eq. (1331) coincide when k is large enough. In 
addition, just like Eq. ([55|. when A/ <C 1, Eq. (1331) will always yield stable solutions, and quasi-static approximation 
will never break down; on the contrary, when A/ is large enough to have a negative coefficient in the first derivative 
term in Eq. (1551) . but not so large to yield hyper-acceleratingly growing solutions, Eq. (1551) will always be a good 
approximation for k ^ kx, regardless the validity of quasi-static approximation. 


C. Recovering Einstein-Hilbert Gravity 

We have shown in Sec. IHI AI that we can reduce the coupled pressure and shear equations into two independent 
second-order differential equations when A/ <C 1/fcy or fc ^ fey. Here, we are going to show that we can even reduce 
these equations further into GR. Under the assumptions. A/ <C l/V. and kXf <C I, all terms proportional to A/, 
including fc^Aj, in Eq. (1551) and Eq. (1251) have only secondary effects compared to TL. Accordingly, both equations will 
simply recover the linear perturbation equation in standard gravity (see e.g. 0 ), 

^ + 3H4'-b (277-f = 0, (36) 

where = $ according to the shear equation, Eq. (El. 

We can easily recognize that Eq. (1361) is independent of scale, k. This differs from the case in f{R) gravity, where 
the evolution of two potentials are scale-dependent. In addition to the scale-independence, Eq. (1561) for standard 
gravity also stands out from its counterparts in f{R) gravity in the following way. According to Eq. (156)) . the impact 
of gravity on the potential only depends implicitly through R. On the contrary, even in the sub-Hubble limit, the 
influence of f{R) gravity is explicitly shown via A/ and its derivatives in Eq. (155)) and Eq. (1551) . In other words, even 
if two f{R) models possess an identical cosmic background evolution, we can still in principle distinguish these two 
models because of their disparity on potentials. This is well known in the previous literature (^ . 


D. Summary of different regimes of linear perturbations in f{R) models 

There is one scale length in f{R) models, A/, which controls the scale on which the modified gravity effects 
contribute. Here, we are going to summarize how this scale length is connected to three different regimes of linear 
perturbations in /(i?) models, and show these connections in Figure El 

(1) GR regime. In the limit Xf ^ I, when f{R) effects vanish, we return to the GR regime, Eq. (1551) . Indeed, 
Aj <C I is a necessary condition to guarantee both kXf and kxXf to be trivial because the minimum value of fcy is 
R, and not negligible. 








kXf 


0 


Reduced second-order f(R) regime / k»kx 
( Eqs. 24, 25, 31-33 ) /' 


k= kv 


kx=:7f 

(Quasi-static regime) 


kx 



Full fourth-order f(R) regime 
(Eqs. 17-20, A1-A4) 


0 


kx A.f 


FIG. 1: Range of applicability of limits of the f(R) equation. In this diagram, both axes increase with distorted scales, and 
we set £ <^ 1. The dashed line represents k = kx, and the dashed-dotted line, kx- When kx\f is not negligible, we enter 
the full fourth-order f{R) regime, and only exact equations are applicable. In the limit fcxA/ <C 1 or fc kx, when k\f is 
not negligible, the reduced second-order f{R) regime is applicable. This regime could be further divided into the quasi-static 
{kx = T-L), and non-quasi-static {kx > 'H) subcategories, which is independent of k, and represented by the dotted line. In the 
limit A/ <C 1, when the regime of modified gravity vanishes, we return to the GR regime. 


(2) Reduced second-order f{R) regime. When kXf is not trivial, but kxXf, compared to I or fcA/, is still negligible, 
we will reach a regime where the exact equations will reduce to decoupled second-order equations, Eqs. (121, (113,®- 
(13311 . which include the so called quasi-Newtonian equation for matter density contrast. Although in the previous 
literature the quasi-static approximation, kx = %, \s assumed in order to derive the quasi-Newtonian equation for 
matter density contrast [H, [3 113, HE HE HE , the quasi-static regime only constitutes part of the reduced second- 
order f{R) regime. Indeed, as long as k ^ kx, no matter the quasi-static approximation breaks down or not, the 
quasi-Newtonian equation for matter density contrast will always be valid. We would also like to emphasize that 
when A/ grows beyond a certain point, the quasi-static approximation will always break down {kx > 'H), and it will 
become more difficult to satisfy the condition, k ^ kx- 

(3) Full fourth-order f{R) regime. When kxXf is not trivial compared to 1 and kXf, it is impossible to further 

reduce the exact equations, Eqs. dH-dlOl), (lATIl-dAl. The only possible way to turn these equations into decoupled 
ones is to recast them into fourth-order differential equations [2J, . We thus call this regime the full fourth-order 

f{R) regime. 


IV. NUMERICAL SOLUTIONS OF LINEAR PERTURBATIONS 

In the previous sections, we have analytically shown that the pressure equation, Eq. dni), and the shear equation, 
Eq. (fT^ . can be used to solve for Newtonian potential, and spatial curvature exactly. We also took these exact 
equations into various limits analytically. In this section, we are going to solve Eq. (jl^^j) and Eq. m numerically 
by turning these two coupled equations into four first-order differential equations, Eq. (TAB-dlll). and apply their 
numerical solutions to check the potential disparity between the exact and quasi-Newtonian solutions at sub-Hubble 
scales, which is suggested by When solving various equations for linear perturbations in f{R) models, we have 
assumed that in the early epoch the evolution of both the Newtonian potential and the Newtonian spatial curvature 
follow that of standard gravity. Accordingly, without losing any generality, we normalize the potentials to ACDM at 
the very beginning, and choose $ = df = I and $ = T = 0ata = 0.001 as our initial conditions. In addition, given a 
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FIG. 2: Example of comparisons of the evolutions of The models, fa{R) = —6flA + , and 

fb{R) = —SHa + are used to compare the evolutions of T solved from different schemes at three scales, fc = 0, 

1, and 600 (in units, c~^'Ho), where Rq = c~^'Hq, and we set c = 1 and Ro = 1- Here, A/o is defined as the scale length. A/, 
at a = 1, and kx ~ R for all k in the both models. These schemes include the exact equations (exact), Eq. (IA1I) - (IA4I) . the 
reduced second-order differential equation (2nd), Eq. (1241) . ACDM, Eq. (1361) . the quasi-Newtonian equation (QN), Eq. (1331) . 
and the super-Hubble equations (SH) in Appendix[Bl Eq. (IB1I) - (IB3I) . For the first model, the solutions of exact, 2nd, and QN 
agree with each other at all scales. Although these solutions also agree with ACDM and SH at fc = 0 and 1, they deviate 
significantly at fc = 600. Also, for this model, ACDM agrees with SH at all scales. For the second model, the solutions of exact, 
2nd, and QN differ at fc = 0 and 1, but still agree at fc = 600. Unlike the first model, ACDM differs from all other solutions, 
even including SH. On the contrary, the SH line is under the exact line at fc = 0. The differences between these two models 
are caused by the fact that /_r and fuR of the second model are 10 times larger than the first one. 


specific f{R) function, we use Eq. ([5]) and Eq. (O to obtain R and R before solving for vF or <&. 


A. Comparisons between solutions of exact equations and various limits 

To solve equations of linear perturbations, we still need to specify f{R) functions, and know the Hubble pa¬ 
rameter beforehand. For this purpose, we apply two specific f{R) models, fa{R) = —Q^A + {R/Ro)~^'^ and 
fb{R) = —5Qa + 10(i?/i?o)~^'^, where Rq = c~^Rq, c = 1, and Rq = 1. We use these two models as exemplars 
because they both yield reasonable background evolution, but at the same time have significant distinguishable linear 
perturbations compared to ACDM. Moreover, since both fn and Jrr of fb{R) are ten times larger than that of fa{R), 
a comparison of linear perturbations between these two models are going to provide us with some hint about how fa 
and fun will affect linear perturbations in f{R) gravity. 

In Fig. [21 we compare the Newtonian potentials solved for the models fa{R) and fb{R) from four exact first-order 
equations, Eq. (Tnt-dlll). from the reduced second-order equations, Eq. (IMl) . and from the quasi-Newtonian second- 
order differential equations, Eq. (l33l) . We plot the solutions at three scales, fc = 0, 1, and 600 (remember k is in units 
of c~^Ro). We can see that all the Newtonian potentials solved from these three different ways concur at fc = 600, 
where k R, and at the early Universe, where A/ <C 1/R. This concurrence continues at the super-Hubble scale, 
fc = 0, as well as the near-Hubble scale, fc = 1, for the model fa{R), but breaks down at the same scales for the model 
fb{R)- As we have shown at the end of Section iHll this difference can be explained by the fact that A/ <C 1/R is still 
valid in the model fa{R) at the recent epoch (A/ = 0.036 at a = 1), but not in the model fb{R) (A/ = 0.179 at a = 1). 

In Section imi we also have shown that under the limit A/ <C 1, linear perturbation in f{R) gravity will recover 
standard gravity. This, however, will break down if the perturbations enter the regime, k~^ \f. This conclusion 

is again supported by Fig. |2j In this figure, the Newtonian solutions solved from the exact equation, the reduced 
second-order equations, and the quasi-Newtonian second-order differential equations all coincide with ACDM for the 
model fa{R) at fc = 0, and 1. This consistency then breaks down at fc = 600 because the effects of f{R) gravity 
caused by fcA/ are no longer negligible at this scale. On the contrary, for the model, fb{R), the effects of /(i?) gravity 
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FIG. 3: Solutions for exact linear perturbations. The exact solutions of and $ are solved from four first-order differential 
equations, Eq. (IAlll - (IA4l) . for the model, fa{R) = —6flA + (R/Ro) between k — 0 and 1000 (in the units, c ^Ro); Ro = 
where we set c = 1, and "Ho = 1. For this f{R) model, 4* and $ at fc = 0 are indistinguishable from their counterparts 
at fc = 1, and ACDM. The combination of the two potentials, -|- $)/2, are shown in the left-lower panel. In the right-lower 
panel, we also plot (’F — 'I>)/2; their counterparts in ACDM are always zero. Because larger k will become compatible with 
the Compton length. A/ earlier than smaller k, f{R) effects will start earlier at sub-Hubble scales. Because effects of f{R) 
are proportional to k^Xj at near or sub-Hubble scales, potentials at smaller scales will deviate from ACDM more dramatically 
than their counterparts at Hubble or super-Hubble scales, where the effects of f{R) gravity are only proportional to A^. 


are nontrivial at all scales because this model possesses the larger value of A/. 

In addition to the approximated equations in the sub-Hubble limit as we have discussed in this paper, offers us 
another approximated equations, Eq. (IB1I) - (IB3I) . for the opposite end of the scale, k = 0. In Fig. [21 we also compare 
the solutions of this approximated equations for k = 0 with our exact solutions solved from Eq. (|M1)-(E1- The 
comparison shows that the Newtonian potentials solved from these two different sets of equations do match at fc = 0, 
but not necessarily in the sub-Hubble limit. All these details prove the consistency between our exact equations and 
Eq. (IB1I) - (IB3I) at fc = 0. 


B. Evolutions of T and $ in f{R) gravity 

In the last section, we compare the solutions of the exact equations with various limits, including A/ -C 1, fc 3> 1, 
kXf ^ 1, and fc = 0. The approximate solutions agreed with the exact solutions in the appropriate limit. In Fig.jSl we 
show that the evolutions ofd), df, (T'-|-<I>)/2 and ('F — <i))/2 solved from the exact equations Eq. jAlll-lAll at several 
scales for the model fa{R)- Unlike the linear perturbation in ACDM, where evolutions of potential are independent of 
fc, both the Newtonian potential and curvature potential at sub-Hubble scales (fc ^ 1) deviate significantly from their 
counterparts at super-Hubble scales. This can be explained by the fact that the effects of f{R) gravity are related to 
kXf at sub-Hubble scales, but are only related to Ay at super-Hubble scales. 

In Fig. O we also show that there exists significant differences between the Newtonian potential and curvature 
potential. At sub-Hubble scales, the Newtonian potentials will mostly grow up from a constant, reach a peak, and 
then decline almost at the same rate for all scales. The curvature potentials at sub-Hubble scales will decline faster 
than their counterparts at super-Hubble scales, oscillate, and eventually, like the Newtonian potentials, decline almost 
at the same rate for all scales. Because for the model fa{R), the exact solutions are identical to the solutions of the 
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a 


a 


FIG. 4: Comparisons of 5 for f{R) = —where Rq = c~^Ro, and we set c = 1, and Ho — 1. The perturbations 
of the matter density in f{R) = —4(i?/i?o)°'®^ are solved from the exact equations (Eq. (1101) and Eq. (IAlll - (IA4ll : exact), the 
quasi-Newtonian equation for 5 (combination of Eq. (1311) and Eq. (1321) ; QN), and ACDM for fc = 0, 1, 10, and, 600 (in units, 
Hole). Although quasi-Newtonian solutions (QN) do not fit the exact solutions (exact) at super-Hubble scale, fc = 0, or around 
Hubble scales, k = 1, there is iro obvious disparity between two solutioirs at near, k — 10, iror sub-Hubble scales, k — 600 Ho/c. 
In this model, kx = H for all k, and the quasi-static assumption has never broken down. 


reduced second-order differential equations, Eq. (l24t and Eq. (l25l) . it is possible to understand why the Newtonian 
and curvature potentials evolve in their particular ways by analyzing Eq. (j24l) and Eq. (12511 themselves. In this way, 
we found that the coefficients of the ih and 4) terms are always positive for the model fa{R), so these terms only 
contribute to damp the evolutions of both potentials. The negative and positive parts of the coefficients of and $ 
contribute to boost and suppress the evolutions of T and d) at sub-Hubble scales. 


C. Issues around quasi-static assumption 

As we have mentioned before, has found a fourth-order differential equation to solve for matter density pertur¬ 
bation exactly. By taking the sub-Hubble limit, they reduce their fourth-order differential equation into a second-order 
one, which they claim is not necessarily identical Eq. (IHTl) combined with Eq. (1321) even at sub-Hubble scales. They 
contribute this disparity to the over-aggressive quasi-static approximation behind the quasi-Newtonian equations 
Eq. (1^ and Eq. (1^ . 

In order to check their claims, we solve for df and $ from the exact equations, Eq. dSl)-® for the model 
fc{R) = —which is used arbitrarily by to show the disparity between the quasi-Newtonian and the 
exact solutions at fc = 600. We then apply these solutions back to Poisson equation, Eq. dini), to obtain matter 
density contrast, 5. We are confident in our solutions because all of them pass the test of consistency by putting 
our solutions back into the shear and pressure equations, Eq. (1181) and Eq. dnj. We compare this 5 with the matter 
density contrast solved from Eq. (IHTl) combined with Eq. (1321) . and plot the results in Eig.|T]for several scales. 

Eig.U] shows that although there exits significant disparity between the exact solutions and the quasi-Newtonian 
ones for fc < 1, this disparity vanishes completely for k > 10, which, of course, includes k = 600. This result conforms 
with our analysis in Section fllll where we show that the difference between the exact solutions and the approximated 
solutions of Eq. Eq. (|25p , cLiici Eq. (|33p circ sol6ly dctcruiincci by th6 competitive terms of cnid k\j iii the 
shear eQuation, Eq. (HU). When k becomes as large as 10 and A/ is not really so large (the maximum value of A/ for 
the model fc{R) is 0.098 between the CMB is formed and today), the reduced second-order differential equations (or 
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quasi-Newtonian equations) will be fairly good approximations. 


V. DISCUSSION AND CONCLUSION 

In this paper, we provide a new insight into how to solve for linearized cosmological perturbations in /(i?) models 
correctly and efficiently. Although in principle, four among two conservation equations and four modified Einstein 
equations of linearized cosmological perturbations are sufficient to solve for the evolutions of $, dl, 0 and S in f{R) 
gravity, in practice it is hard to solve these coupled equations simultaneously. Many studies circumvent this problem by 
tak ing insi ghtf ul assumptions for some specific circumstance, such as quasi-static or sub-Hubble approximations |18l . 
flfll . [U, [25l427j |. In this paper, we analytically prove that without losing any information, the pressure equation, 
Eq. (UHl), and the shear equation, Eq. (fT^ . alone are sufficient to solve for $ and d', which although is obvious in GR, 
has never been shown explicitly in the f{R) literature. After obtaining evolutions of dl and $, we can easily derive <5 
and 9 from the Poisson equation, Eq. (HOD, and the momentum equation, Eq. (USD- 

One of the advantages to apply the pressure equation, Eq. (HID, and the shear equation, Eq. dUD, rather than the 
equations used in (20l - l^ 1^ , is that we can easily understand how Eq. (1151) and Eq. (1151) will reduce into the decoupled 
second-order differential equations, Eq. (l24l) and Eq. (USD, which themselves can be applied for the consistent analysis 
of the quasi-static approximation. In our reductions, we naturally find a Compton length. Ay, to characterize the 
effects of f{R) gravity in linear perturbations. We conclude how the exact equations, Eq. (fT51) - (lT5D . fit into the 
reduced second-order differential equations Eq. (1M1) - (1551) . quasi-Newtonian equations, Eq. (1551) - (1551) . and Einstein 
gravity, Eq. ([36D, solely depends on the effects of the A/ and kXj terms in the shear equation, Eq. dSD. We reach the 
same conclusion of [2,1^ from a different approach. 

We also find that A / plays an important role in the quasi-stati c appro ximation: when Ay is small, which is applicable 
for most if not all observational viable/(i?) models lO Il4. fla. lUU - H^ . quasi-static approximation will always be 
valid. Our Ending agrees well with |3l| . where the quasi-static approximation has been proved valid for |//j| = 
10“"^ — 10“®. Because small Ay will generically result in ACDM-like background evolutions, our conclusion is also 
consistent with [23, which find that ACDM-like backgrounds will guarantee the applicability of quasi-static 
approximation. However, we would like to emphasize that although a small Compton length will guarantee a proximity 
of background evolutions to ACDM, it is not always correct the other way around. As shown in [2^, different f{R) 
models might lead to the same background evolution that is identical to ACDM. Our analysis shows that it is Ay 
rather than a proximity of background evolutions to ACDM that will guarantee the quasi-static approximation. This 
finding might even have a deep connection to the conclusion of [30l| . which shows that the quasi-static behaviors are 
related to the sound horizon of dark energy or modified gravity; that is, the nature of models themselves. We would 
like to explore this plausible connection in the future. 

In addition to quasi-static approximation, we also discuss the role of sub-Hubble limit in taking approximations 
of the exact equations, Eq. (fT8l) - (IT9l) . Our analysis also shows that whether or not quasi-static approximation beaks 
down, as long as Ay is not extremely large we can always find a k that is large enough to guarantee the reduced 
second-order differential equations, Eq. dMD-dlSD, and quasi-Newtonuan equations, Eq. (1551) - (1551) . Indeed, even for 
the model that found problematic in their quasi-Newtonian solutions at fc = 600, we cannot find a ny obvious 
disparity between the exact and quasi-Newtonian solutions up to fc = 10. Our conclusion contradicts |24| . but is 
consistent with [32|. 

Generally speaking, the pressure equation, Eq. dUD, and the shear equation, Eq. dUD are only based on the 
assumption of a flat RW metric and a matter-dominated era. Therefore combined with Poisson equation, these two 
equations shall be legitimate to be applied to inve stig ate oscillating solutions of matter density perturbations. Since 
our analysis casts doubts in some conclusions in |24| . it will be interesting to re-analyze these oscillating solutions 
from our approach, and compare them with the results in (^ . We leave the exploration of this study to the near 
future. 
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Appendix A: Exact Decoupled Equations for Linear Perturbations in f{R) Gravity 


We have seen the pressure equation, Eq. dn, and the shear equation, Eq. (IT^ . are sufficient to solve for <i> and dl 
exactly. Practically speaking, in order to solve for these two potentials, we might apply numerical method by turning 
Eq. (fTSl) and Eq. (fTOl) into four coupled first-order differential equations: 


where. 


^ = 0, (Al) 

$ = fl, (A2) 

0 = 0([),0,d/,$), (A3) 

f) = Jo((5,0,vI/,$), (A4) 


© = 


' aHl + fn) ^ 2^2 _ ^ ^ nfn + /it 

" (ijRR 3 1 + 




= - 


' aHl + fR) 
^/rr 

^ 0^(1 + Jr) 
^/rr 


\e-n 




■)*+( 


^'HJr + 2fR - a^Si rCp 
1 + /it 

^-il±M + h^-2in + n^ 

OjRR 3 



(AS) 


(A6) 


Although in practical Eq. (EH) -Eq. (IA6I) are applied to solve for the Newtonian potential, di, and the curvature 
potential, $, numerically, the pressure equation, Eq. m, and the shear equation, Eq. di), themselves will provide 
us with a more intuitive sense about how the effects of f{R) gravity will contribute to the linearized cosmological 
perturbations. 


Appendix B: The Super-Hubble limit 


In their paper about the large scale structure of f{R) models, 
at fc = 0: 


offers us a simple relation between 'k and $ 


$ -I- -I- 

B-1 ~ n{B-i) ■ 


(Bl) 


They also provide us an equation that can exactly describe the linear perturbation at the super-Hubble scale, k = 0. 
In our notation this equation appears as 


where 


- / , H 2 nn-2nn'^ b b Bn\- 

\ n a nn^-2'H^ an{l-B) a an^J 

/-n n iin - 2nn^ - b ^ 

+ (—+ S7+ HW-2W +;(T3b)) 


JrR RR _ IrR p R 
1 + fR aR' -R~ 1 + fR R-R'^ 


(B2) 


(B3) 


is a specific parameter used by [l^, and B^/"^ is called Compton length by 0- Analytically, there seems no obvious 
way to prove that the pressure and shear equations can be reduced into Eq. (IBlIl and Eq. (IB2I1 at fe = 0. However, our 
numerical solutions in Sec. IIVI show that the Newtonian potentials from both sides do concur at this very super-Hubble 
scale. 
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